clear; close all; clc;

%% Load GDP data

GDPdata = xlsread('Data_MATLAB.xlsx','GDP_PC','D1:BS15');

% Time series of the log of GDP per capita for each country
year =     GDPdata(1,:)'; 
ARG  = log(GDPdata(2,:))';
BRA  = log(GDPdata(3,:))';
ITA  = log(GDPdata(4,:))';
POR  = log(GDPdata(5,:))';
SPA  = log(GDPdata(6,:))';
% One-year differences
dyARG  = ARG(2:end) - ARG(1:end-1);
dyBRA  = BRA(2:end) - BRA(1:end-1);
dyITA  = ITA(2:end) - ITA(1:end-1);
dyPOR  = POR(2:end) - POR(1:end-1);
dySPA  = SPA(2:end) - SPA(1:end-1);
dy     = [dyITA,dyPOR,dySPA,dyARG,dyBRA];
year   = year(2:end);
ncountries = 5;

%% Set-ups for estimation

% Output: tables with mean, 5th and 95th percentiles of each parameter for
% each country
Ntheta     = 5;                             % number of parameters
table_mean = zeros(ncountries,Ntheta);      % means
table_p05  = zeros(ncountries,Ntheta);      % 5th percentiles
table_p95  = zeros(ncountries,Ntheta);      % 95th percentiles
NS = 2;  % number of states for Markov process
Nx = 2;  % length of unobserved vector
NC = 1;  % number of countries in each estimation (we estimate each country separately)
cnames = {'Italy','Portugal','Spain','Argentina','Brazil'}; 
% The scale parameters target rejection rates between 60% and 70%
scale_parameter = zeros(ncountries,1);
scale_parameter(1) = 0.015;      % Italy
scale_parameter(2) = 0.038;      % Portugal
scale_parameter(3) = 0.012;      % Spain
scale_parameter(4) = 0.030;      % Argentina
scale_parameter(5) = 0.006;      % Brazil
MH.CH    = 1;                    % number of chains to simulate
MH.Nburn = 25000;                % number of iterations for burn-in sample
MH.Nsim  = 225000;               % number of iterations to simulate draws
Nburn2 = 5000;                   % number of iterations to burn when computing the moments of the posterior distribution
MH.NB    = MH.Nburn + MH.Nsim;   % number of total iterations
MH.ext   = 10;                   % save a draw every "ext" elements of the chain
MH.print = 1;                    % print results if = 1 
MH.show  = 1000;
MH.save  = 1000;                

% Defining prior distributions - uniform
prior.p11lb = 0.05;    prior.p11ub = 0.99;  % bounds on the prior distribution for p11 (persistence of gL)
prior.p22lb = 0.05;    prior.p22ub = 0.99;  % bounds on the prior distribution for p22 (persistence of gH)
prior.mu1lb = -0.1;    prior.mu1ub = 0.1;   % bounds on the prior distribution for mu1 (gL)
prior.Dmulb =  0.0;    prior.Dmuub = 0.1;   % bounds on the prior distribution for m2-mu1; (gH = gL + (mu2-mu1))
prior.sglb  =  0.001;  prior.sgub  = 0.075; % bounds on the prior distribution for sg (standard deviation of transitory shock)

%% Loop over countries

for iic = 1:ncountries

if iic==1
     rng(2) % we got a corner solution with rng(1) for the case of Italy. Estimates were not converging.
else
    rng(1)
end

ydata = dy(:,iic);
cname  = cnames(iic);

% We select the initial point from the solution of the maximum likelihood estimation (MLE)
mu100 = mean(ydata(ydata<mean(ydata)));  % initial guess for gL to compute the MLE
mu200 = mean(ydata(ydata>=mean(ydata))); % initial guess for gH to compute the MLE
Dmu00 = mu200 - mu100;                   % initial guess for (gH-gL)
esim = ydata - (ydata<mean(ydata))*mu100 - (ydata>=mean(ydata))*mu200;
sg00 = std(esim)/2;                      % initial guess for sg (standard deviation of transitory shock) to compute the MLE
N11   = sum( (ydata(1:end-1)<mean(ydata)).*(ydata(2:end)<mean(ydata)) );
N1    = sum( (ydata<mean(ydata)) );
p1100 = N11/N1;                          % initial guess for p11 (persistence of gL) to compute the MLE
N22   = sum( (ydata(1:end-1)>=mean(ydata)).*(ydata(2:end)>=mean(ydata)) );
N2    = sum( (ydata>=mean(ydata)) );
p2200 = N22/N2;                          % initial guess for p22 (persistence of gH) to compute the MLE
% in the maximization/minimization we choose the weight between the lower bounds and upper bounds in the prior distribution of each parameter
% the weights must be between 0 and 1, so we choose log(weight/(1-weight));
% vector of weights in our initial guess:
htheta = [ (p1100-prior.p11ub)/(prior.p11lb-prior.p11ub); ...
           (p2200-prior.p22ub)/(prior.p22lb-prior.p22ub); ...
           (mu100-prior.mu1ub)/(prior.mu1lb-prior.mu1ub); ...
           (Dmu00-prior.Dmuub)/(prior.Dmulb-prior.Dmuub); ...
           (sg00 -prior.sgub )/(prior.sglb -prior.sgub )]; 
% initial guess for log(weight/(1-weight))
xx00    = log(htheta./(1-htheta));      
% Maximum Likelihood Estimation (MLE)
tr = 1; setmin = 1;
%options = optimset('Display','iter');
[xx_out,fval] = fminsearch(@(xx)posterior_fun_MA(xx,ydata,prior,NS,Nx,NC,tr,setmin),xx00);
% xx_out is our initial guess for log(weights/(1-weights)) for each
% parameter

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computing the MH Algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xx0      = xx_out;              % initial guess from the MLE
scale = scale_parameter(iic);   % scale parameter targeting rejection rate between 60% and 70%. 
% Covariance Matrix for draws --> update after burn-in periods:
SIGMA_xx    = scale*diag(max(abs(xx0),0.1));    
SIGMA_xx    = 0.5*(SIGMA_xx + SIGMA_xx');
xx_chain    = 189*ones(Ntheta,MH.NB);
lnL_chain   = 189*ones(MH.NB,1);
RAT_store   = 189*ones(MH.NB,1);
uu_draw     = unifrnd(0,1,MH.NB,1);  % for the acceptance/rejection decision
tr          = 1;   % transformation for posterior fun
setmin      = 0;   % = 1 if -PS in posterior fun  

% Initialize Chain
ib = 1;
xx_chain(:,ib) = xx0;
lnLch          = posterior_fun_MA(xx0,ydata,prior,NS,Nx,NC,tr,setmin);
lnL_chain(ib)  = lnLch;

% Markov Chain simulation
acc = 0; rej = 0; ibtot = 0; ss = MH.save; show = MH.show;
ch = 1; 
for ib = 2:MH.NB
    
    ibtot = ibtot + 1;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%    New draw Chain   %%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%   
    lnLj      = lnL_chain(ib-1);
    xx_new    = mvnrnd(xx_chain(:,ib-1),SIGMA_xx);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%   Update the Chain   %%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    lnL_new       = posterior_fun_MA(xx_new,ydata,prior,NS,Nx,NC,tr,setmin);                
    RAT_new       = exp(lnL_new-lnLj);         % exp(lnL_new)/exp(lnLj);        
    RAT_store(ib) = RAT_new;    
    if RAT_new > uu_draw(ib)  % this is a good theta draw -->  update
        xx_chain(:,ib) = xx_new;
        lnL_chain(ib)  = lnL_new;        
        acc = acc + 1;
    else
        xx_chain(:,ib) = xx_chain(:,ib-1);
        lnL_chain(ib)  = lnL_chain(ib-1);        
        rej = rej + 1;
    end
    
    if MH.print==1 && ch == 1   
        if show == 0
        disp(['MH   ',cell2mat(cnames(iic)),'  ib = ', num2str(ib),', acc = ', num2str(acc/(acc+rej),3),', rej = ', num2str(rej/(acc+rej),3), ', lnL = ', num2str(lnL_chain(ib),3),', eff = ', num2str((acc+rej)/ibtot,3)])
        show = MH.show;
        end
    end
    
    if ib == MH.Nburn  % reset acc and rej count
        disp('************************************************')
        disp('*****  resetting acc, rej and ibtot count  *****')
        disp('************************************************')
        acc = 0; rej = 0; ibtot = 0;
    end
    
    show = show - 1;
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Posterior computations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% transform xx into theta
theta_chain = xx_chain;
alphap11 = exp(xx_chain(1,:))./(1+exp(xx_chain(1,:))); theta_chain(1,:)   = alphap11.*prior.p11lb + (1-alphap11).*prior.p11ub; % p1
alphap22 = exp(xx_chain(2,:))./(1+exp(xx_chain(2,:))); theta_chain(2,:)   = alphap22.*prior.p22lb + (1-alphap22).*prior.p22ub; % p2
alphamu1 = exp(xx_chain(3,:))./(1+exp(xx_chain(3,:))); theta_chain(3,:)   = alphamu1.*prior.mu1lb + (1-alphamu1).*prior.mu1ub; % gH
alphaDmu = exp(xx_chain(4,:))./(1+exp(xx_chain(4,:))); theta_chain(4,:)   = alphaDmu.*prior.Dmulb + (1-alphaDmu).*prior.Dmuub; % gH-gL
alphasg  = exp(xx_chain(5,:))./(1+exp(xx_chain(5,:))); theta_chain(5,:)   = alphasg .*prior.sglb  + (1-alphasg) .*prior.sgub;  % sg
Dmu_chain        = theta_chain(4,:);           
theta_chain(4,:) = theta_chain(3,:) + Dmu_chain;  % transforming (gH-gL) into gH
% save entire chain
theta_chain_lng = theta_chain;
% discard MH.Nburn initial draws, and then keep one every MH.ext draws
in          = MH.Nburn+1:MH.ext:MH.Nsim;
theta_chain = theta_chain_lng(:,in);

% Computing mean and 5th and 95th percentiles
table_mean(iic,:) = mean(theta_chain(:,Nburn2:end),2)';
table_p05(iic,:)  = prctile(theta_chain(:,Nburn2:end)',5);
table_p95(iic,:)  = prctile(theta_chain(:,Nburn2:end)',95);

end


filename = 'Table_1.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)

% Table 1
fprintf('------------------------------------------------------------- \n')
fprintf('TABLE 1: Posterior distributions \n')
fprintf('------------------------------------------------------------- \n')
fprintf('           \t   gL   \t  gH   \t   pL \t  pH   \t  sg   \n')
for iic=1:ncountries
fprintf([cell2mat(cnames(iic)),'\n'])
fprintf(' mean     \t  %5.3f \t %5.3f \t %5.3f \t %5.3f \t %5.3f \n', [table_mean(iic,3),table_mean(iic,4),table_mean(iic,1),table_mean(iic,2),table_mean(iic,5)])
fprintf(' p05     \t  %5.3f \t %5.3f \t %5.3f \t %5.3f \t %5.3f \n', [table_p05(iic,3),table_p05(iic,4),table_p05(iic,1),table_p05(iic,2),table_p05(iic,5)])
fprintf(' p95     \t  %5.3f \t %5.3f \t %5.3f \t %5.3f \t %5.3f \n', [table_p95(iic,3),table_p95(iic,4),table_p95(iic,1),table_p95(iic,2),table_p95(iic,5)])
fprintf('\n')
end
fprintf('------------------------------------------------------------- \n')

diary off
